// adapted by pbm from make_meson.cpp by mdp, Aug. 05
// #define PARALLEL
#include "fermiqcd.h"
#include "fermiqcd_ffts.h"
#include "math.h"
#include "fermiqcd_fermilab_action.h"

//#define USE_S

int main(int argc, char **argv) {
  mpi.open_wormholes(argc,argv);
  define_base_matrices("FERMILAB");
  
  // Setting parameters here:
  coefficients gauge;
  gauge["beta"]=6.0;

  coefficients light_quark;
  light_quark["kappa_s"]=0.1;
  light_quark["kappa_t"]=0.1;
  light_quark["c_{sw}"]=0;
  light_quark["c_E"]=0;
  light_quark["c_B"]=0;

  light_quark["r_s"]=0;
  light_quark["r_t"]=0;

  light_quark["mac_3"]=-1./6.;
  light_quark["alpha_2"]=-1./6.;
  light_quark["alpha_4"]=-1./4.;

  double absolute_precision=1e-23;
  double relative_precision=1e-12;
  
  /*
#ifdef SSE2
default_fermi_action=FermiCloverActionSSE2::mul_Q;
#else
  default_fermi_action=FermiCloverActionFast::mul_Q;
#endif
  */

  default_fermi_action=FermiFermilabActionNew::mul_Q;
  
  default_fermi_inverter=MinimumResidueInverter<fermi_field,gauge_field>;

  int nt=16, ns=6, nc=3, t,a,b,i,j,n;
  double alpha=0.0;
  double pi=asin(1.)*2.;

  int mybox[]={nt,ns,ns,ns};

  // creating the fields
  generic_lattice mylattice(4,mybox);
  gauge_field U(mylattice,nc);
  fermi_field psi(mylattice,nc);
  fermi_field phi(mylattice,nc);
  fermi_field phi_ft(mylattice,nc);
  site x(mylattice),x0(mylattice),x1(mylattice);
  mdp_array<Complex,1> C2(nt);

  // creating one gauge configuration

  set_cold(U);
  // WilsonGaugeAction::heatbath(U,gauge,20); // 2 iterations
  compute_em_field(U); // THIS MISSINS CAUSES RUNTIME ERROR

  int ms=1;
  int mc=1;

#ifdef USE_S
  
  fermi_propagator S(mylattice,nc);

  generate(S,U,light_quark,absolute_precision,relative_precision);
  
  for(t=0; t<nt; t++) C2(t)=0;

  forallsites(x)
    for(a=0; a<ms; a++)
      for(b=0; b<ms; b++)
	for(i=0; i<mc; i++)
	  for(j=0; j<mc; j++){
	    C2(x(0))+=S(x,a,b,i,j)  ;   // *conj(S(x,a,b,i,j));
	  }

#else

  for(t=0; t<nt; t++) C2(t)=0.;

  for(a=0; a<ms; a++)
    for(i=0; i<mc; i++) {      
      psi=0;
      if(on_which_process(mylattice,0,0,0,0)==ME) {
	x.set(0,0,0,0);
	psi(x,a,i)=1;
      }
      psi.update();

      // make light prop phi
      mul_invQ(phi,psi,U,light_quark,absolute_precision,relative_precision); 

      fermi_field_fft(phi_ft, phi, +1);
      phi_ft*=pow((1.*ns),1.5);

//       forallsites(x)
// 	for(b=0; b<ms; b++)
// 	  for(j=0; j<mc; j++)
// 	    C2(x(0))+=phi(x,b,j) *cos(x(1)*2.*pi/ns);

    }
  
#endif

  mpi.add(C2.address(),nt);
  for(t=0; t<nt; t++) {
    if(on_which_process(mylattice,t)==ME) { 
      x0.set(t,0,0,0);
      x1.set(t,1,0,0);
      printf("%i\t%e\t%e\n", t, real(phi_ft(x0,0,0)), real(phi_ft(x1,0,0)));
      //     printf("%i\t%e\t%e\n", t, real(C2(t)), imag(C2(t)));
    }
    mpi.barrier(); // to make sure you print sites in order in parallel
  }

  mpi.close_wormholes();
};

